{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "extraordinary-scenario",
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx\n",
    "import scipy.sparse as sprs\n",
    "import numpy as np\n",
    "#import scipy.sparse\n",
    "import scipy.sparse as sprs\n",
    "import copy \n",
    "import random as r_n\n",
    "import pylab\n",
    "from random import randint \n",
    "from matplotlib import mlab\n",
    "from scipy.stats import bernoulli\n",
    "from scipy.optimize import minimize\n",
    "import pandas as pd\n",
    "from scipy.sparse import csr_matrix, lil_matrix, hstack, vstack\n",
    "import random\n",
    "from tqdm.auto import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import preprocessing \n",
    "from sklearn.preprocessing import normalize\n",
    "from random import choice\n",
    "from scipy.stats import binom\n",
    "import os\n",
    "from sklearn.model_selection import ParameterGrid\n",
    "from multiprocessing import Process\n",
    "from IPython.display import clear_output\n",
    "from statsmodels.stats.proportion import proportion_confint\n",
    "from IPython.display import clear_output\n",
    "import math\n",
    "from scipy import integrate\n",
    "import seaborn as sns\n",
    "from IPython.display import clear_output\n",
    "import time\n",
    "from sympy import Symbol\n",
    "from sympy import *\n",
    "from scipy.optimize import minimize, rosen, rosen_der"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interior-ancient",
   "metadata": {},
   "source": [
    "# Diversity of opinions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "threatened-parks",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ops1</th>\n",
       "      <th>ops2</th>\n",
       "      <th>ops3</th>\n",
       "      <th>fr_ops1_avg</th>\n",
       "      <th>fr_ops1_std</th>\n",
       "      <th>fr_ops2_avg</th>\n",
       "      <th>fr_ops2_std</th>\n",
       "      <th>fr_ops3_avg</th>\n",
       "      <th>fr_ops3_std</th>\n",
       "      <th>fr_number</th>\n",
       "      <th>SL-Friends</th>\n",
       "      <th>L-Friends</th>\n",
       "      <th>M-Friends</th>\n",
       "      <th>C-Friends</th>\n",
       "      <th>SC-Friends</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.514611</td>\n",
       "      <td>0.514624</td>\n",
       "      <td>0.512251</td>\n",
       "      <td>0.410485</td>\n",
       "      <td>0.261048</td>\n",
       "      <td>0.467247</td>\n",
       "      <td>0.251042</td>\n",
       "      <td>0.457429</td>\n",
       "      <td>0.252795</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.186235</td>\n",
       "      <td>0.186235</td>\n",
       "      <td>0.186192</td>\n",
       "      <td>0.544541</td>\n",
       "      <td>0.227728</td>\n",
       "      <td>0.576840</td>\n",
       "      <td>0.238795</td>\n",
       "      <td>0.577952</td>\n",
       "      <td>0.239873</td>\n",
       "      <td>7.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.604282</td>\n",
       "      <td>0.603247</td>\n",
       "      <td>0.603549</td>\n",
       "      <td>0.333697</td>\n",
       "      <td>0.167857</td>\n",
       "      <td>0.308934</td>\n",
       "      <td>0.144153</td>\n",
       "      <td>0.286581</td>\n",
       "      <td>0.132653</td>\n",
       "      <td>6.0</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.608380</td>\n",
       "      <td>0.614673</td>\n",
       "      <td>0.637561</td>\n",
       "      <td>0.410863</td>\n",
       "      <td>0.269629</td>\n",
       "      <td>0.413971</td>\n",
       "      <td>0.264415</td>\n",
       "      <td>0.413348</td>\n",
       "      <td>0.264016</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.729980</td>\n",
       "      <td>0.722895</td>\n",
       "      <td>0.726959</td>\n",
       "      <td>0.605165</td>\n",
       "      <td>0.108791</td>\n",
       "      <td>0.616444</td>\n",
       "      <td>0.104925</td>\n",
       "      <td>0.614579</td>\n",
       "      <td>0.106953</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       ops1      ops2      ops3  fr_ops1_avg  fr_ops1_std  fr_ops2_avg  \\\n",
       "0  0.514611  0.514624  0.512251     0.410485     0.261048     0.467247   \n",
       "1  0.186235  0.186235  0.186192     0.544541     0.227728     0.576840   \n",
       "2  0.604282  0.603247  0.603549     0.333697     0.167857     0.308934   \n",
       "3  0.608380  0.614673  0.637561     0.410863     0.269629     0.413971   \n",
       "4  0.729980  0.722895  0.726959     0.605165     0.108791     0.616444   \n",
       "\n",
       "   fr_ops2_std  fr_ops3_avg  fr_ops3_std  fr_number  SL-Friends  L-Friends  \\\n",
       "0     0.251042     0.457429     0.252795        7.0           2          0   \n",
       "1     0.238795     0.577952     0.239873        7.0           1          0   \n",
       "2     0.144153     0.286581     0.132653        6.0           2          2   \n",
       "3     0.264415     0.413348     0.264016        3.0           1          0   \n",
       "4     0.104925     0.614579     0.106953        3.0           0          0   \n",
       "\n",
       "   M-Friends  C-Friends  SC-Friends  \n",
       "0          3          1           1  \n",
       "1          1          5           0  \n",
       "2          2          0           0  \n",
       "3          1          1           0  \n",
       "4          2          1           0  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv('main_data.csv')\n",
    "n = data.shape[0]\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "confirmed-addition",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "113d9af194f84e5eb5497b8b1ededc79",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=10.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "42465bd0f47144b992edd25c619f1e2f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=10.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "ths_op = np.arange(0, 1.1, 0.1)  #m=10\n",
    "#ths_op = np.array([0, 0.33, 0.66, 1])  #m=3\n",
    "#ths_op = np.array([0, 0.5, 1])  #m=2\n",
    "\n",
    "\n",
    "#If use the first opinion snapshot\n",
    "op = 'ops1'\n",
    "fr_op = 'fr_ops1_avg'\n",
    "\n",
    "\n",
    "#If use the second opinion snapshot\n",
    "#op = 'ops2'\n",
    "#fr_op = 'fr_ops2_avg'\n",
    "\n",
    "\n",
    "#If use the third opinion snapshot\n",
    "#op = 'ops3'\n",
    "#fr_op = 'fr_ops3_avg'\n",
    "\n",
    "\n",
    "masks = []\n",
    "masks_fr = []\n",
    "\n",
    "for i in tqdm(range(len(ths_op) - 1)):\n",
    "    \n",
    "    if i < len(ths_op)  - 2:\n",
    "        mask_temp_curr_pos = (ths_op[i] <= data[op]) & (data[op] < ths_op[i+1])\n",
    "        mask_temp_fr_pos = (ths_op[i] <= data[fr_op]) & (data[fr_op] < ths_op[i+1])\n",
    "    else:\n",
    "        mask_temp_curr_pos = (ths_op[i] <= data[op]) & (data[op] <= ths_op[i+1])\n",
    "        mask_temp_fr_pos = (ths_op[i] <= data[fr_op]) & (data[fr_op] <= ths_op[i+1])\n",
    "   \n",
    "    masks.append(mask_temp_curr_pos)\n",
    "    masks_fr.append(mask_temp_fr_pos)\n",
    "    \n",
    "    \n",
    "populations = np.zeros((len(masks), len(masks_fr)))\n",
    "\n",
    "for i in tqdm(range(len(masks))):\n",
    "    for j in range(len(masks_fr)):\n",
    "        \n",
    "        populations[i,j] = sum(masks[i] & masks_fr[j])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "genuine-municipality",
   "metadata": {},
   "source": [
    "# Table A1 in Appendix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "circular-midnight",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>$x_1$</th>\n",
       "      <th>$x_2$</th>\n",
       "      <th>$x_3$</th>\n",
       "      <th>$x_4$</th>\n",
       "      <th>$x_5$</th>\n",
       "      <th>$x_6$</th>\n",
       "      <th>$x_7$</th>\n",
       "      <th>$x_8$</th>\n",
       "      <th>$x_9$</th>\n",
       "      <th>$x_10$</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>$x_1$</th>\n",
       "      <td>371.0</td>\n",
       "      <td>699.0</td>\n",
       "      <td>3833.0</td>\n",
       "      <td>21084.0</td>\n",
       "      <td>24955.0</td>\n",
       "      <td>6934.0</td>\n",
       "      <td>1332.0</td>\n",
       "      <td>376.0</td>\n",
       "      <td>147.0</td>\n",
       "      <td>44.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_2$</th>\n",
       "      <td>246.0</td>\n",
       "      <td>486.0</td>\n",
       "      <td>2761.0</td>\n",
       "      <td>19626.0</td>\n",
       "      <td>27570.0</td>\n",
       "      <td>7475.0</td>\n",
       "      <td>1359.0</td>\n",
       "      <td>409.0</td>\n",
       "      <td>124.0</td>\n",
       "      <td>45.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_3$</th>\n",
       "      <td>323.0</td>\n",
       "      <td>595.0</td>\n",
       "      <td>3307.0</td>\n",
       "      <td>27307.0</td>\n",
       "      <td>49928.0</td>\n",
       "      <td>15408.0</td>\n",
       "      <td>2750.0</td>\n",
       "      <td>745.0</td>\n",
       "      <td>217.0</td>\n",
       "      <td>75.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_4$</th>\n",
       "      <td>624.0</td>\n",
       "      <td>1002.0</td>\n",
       "      <td>5064.0</td>\n",
       "      <td>45533.0</td>\n",
       "      <td>105023.0</td>\n",
       "      <td>36341.0</td>\n",
       "      <td>6148.0</td>\n",
       "      <td>1546.0</td>\n",
       "      <td>544.0</td>\n",
       "      <td>151.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_5$</th>\n",
       "      <td>1318.0</td>\n",
       "      <td>2111.0</td>\n",
       "      <td>9915.0</td>\n",
       "      <td>90305.0</td>\n",
       "      <td>277315.0</td>\n",
       "      <td>107562.0</td>\n",
       "      <td>18666.0</td>\n",
       "      <td>4586.0</td>\n",
       "      <td>1575.0</td>\n",
       "      <td>405.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_6$</th>\n",
       "      <td>825.0</td>\n",
       "      <td>1253.0</td>\n",
       "      <td>4932.0</td>\n",
       "      <td>41661.0</td>\n",
       "      <td>193936.0</td>\n",
       "      <td>100509.0</td>\n",
       "      <td>20220.0</td>\n",
       "      <td>5237.0</td>\n",
       "      <td>1833.0</td>\n",
       "      <td>439.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_7$</th>\n",
       "      <td>470.0</td>\n",
       "      <td>614.0</td>\n",
       "      <td>2142.0</td>\n",
       "      <td>14414.0</td>\n",
       "      <td>78395.0</td>\n",
       "      <td>56749.0</td>\n",
       "      <td>14657.0</td>\n",
       "      <td>4079.0</td>\n",
       "      <td>1293.0</td>\n",
       "      <td>350.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_8$</th>\n",
       "      <td>245.0</td>\n",
       "      <td>336.0</td>\n",
       "      <td>1206.0</td>\n",
       "      <td>6922.0</td>\n",
       "      <td>37826.0</td>\n",
       "      <td>34854.0</td>\n",
       "      <td>11471.0</td>\n",
       "      <td>3495.0</td>\n",
       "      <td>1170.0</td>\n",
       "      <td>319.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_9$</th>\n",
       "      <td>179.0</td>\n",
       "      <td>222.0</td>\n",
       "      <td>722.0</td>\n",
       "      <td>3451.0</td>\n",
       "      <td>17376.0</td>\n",
       "      <td>19377.0</td>\n",
       "      <td>8136.0</td>\n",
       "      <td>2703.0</td>\n",
       "      <td>1006.0</td>\n",
       "      <td>273.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$x_10$</th>\n",
       "      <td>58.0</td>\n",
       "      <td>63.0</td>\n",
       "      <td>254.0</td>\n",
       "      <td>1305.0</td>\n",
       "      <td>5553.0</td>\n",
       "      <td>5716.0</td>\n",
       "      <td>2679.0</td>\n",
       "      <td>1063.0</td>\n",
       "      <td>457.0</td>\n",
       "      <td>119.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         $x_1$   $x_2$   $x_3$    $x_4$     $x_5$     $x_6$    $x_7$   $x_8$  \\\n",
       "$x_1$    371.0   699.0  3833.0  21084.0   24955.0    6934.0   1332.0   376.0   \n",
       "$x_2$    246.0   486.0  2761.0  19626.0   27570.0    7475.0   1359.0   409.0   \n",
       "$x_3$    323.0   595.0  3307.0  27307.0   49928.0   15408.0   2750.0   745.0   \n",
       "$x_4$    624.0  1002.0  5064.0  45533.0  105023.0   36341.0   6148.0  1546.0   \n",
       "$x_5$   1318.0  2111.0  9915.0  90305.0  277315.0  107562.0  18666.0  4586.0   \n",
       "$x_6$    825.0  1253.0  4932.0  41661.0  193936.0  100509.0  20220.0  5237.0   \n",
       "$x_7$    470.0   614.0  2142.0  14414.0   78395.0   56749.0  14657.0  4079.0   \n",
       "$x_8$    245.0   336.0  1206.0   6922.0   37826.0   34854.0  11471.0  3495.0   \n",
       "$x_9$    179.0   222.0   722.0   3451.0   17376.0   19377.0   8136.0  2703.0   \n",
       "$x_10$    58.0    63.0   254.0   1305.0    5553.0    5716.0   2679.0  1063.0   \n",
       "\n",
       "         $x_9$  $x_10$  \n",
       "$x_1$    147.0    44.0  \n",
       "$x_2$    124.0    45.0  \n",
       "$x_3$    217.0    75.0  \n",
       "$x_4$    544.0   151.0  \n",
       "$x_5$   1575.0   405.0  \n",
       "$x_6$   1833.0   439.0  \n",
       "$x_7$   1293.0   350.0  \n",
       "$x_8$   1170.0   319.0  \n",
       "$x_9$   1006.0   273.0  \n",
       "$x_10$   457.0   119.0  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.set_printoptions(suppress=True)\n",
    "pd.DataFrame(populations, \n",
    "             columns = [f'$x_{i}$' for i in range(1, len(ths_op) - 1+1)], \n",
    "             index = [f'$x_{i}$' for i in range(1, len(ths_op) - 1+1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ancient-wagon",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(suppress=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "nominated-projector",
   "metadata": {},
   "source": [
    "# Estimation of transition matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "latter-mongolia",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('main_data.csv')\n",
    "n = data.shape[0]\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "liquid-plaza",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ths_op = np.arange(0, 1.1, 0.1)  #m=10\n",
    "#ths_op = np.array([0, 0.33, 0.66, 1])  #m=3\n",
    "ths_op = np.array([0, 0.5, 1])  #m=2\n",
    "\n",
    "\n",
    "#If use first two opinion snapshots (first iteration of opinion change)\n",
    "op1 = 'ops1'\n",
    "op2 = 'ops2'\n",
    "fr_op = 'fr_ops1_avg'\n",
    "\n",
    "\n",
    "#If use second and third opinion snapshots (second iteration of opinion change)\n",
    "#op1 = 'ops2'\n",
    "#op2 = 'ops3'\n",
    "#fr_op = 'fr_ops2_avg'\n",
    "\n",
    "\n",
    "\n",
    "op_change = data[op2] - data[op1]\n",
    "anchor = data[fr_op] - data[op1]\n",
    "\n",
    "mask_move = abs(op_change) > 0.05\n",
    "\n",
    "mask_neg = op_change * anchor < 0\n",
    "\n",
    "mask_pos = op_change * anchor > 0\n",
    "mask_skip = mask_pos & (abs(op_change) > abs(anchor))\n",
    "mask_no_skip = mask_pos & (abs(op_change) <= abs(anchor))\n",
    "\n",
    "masks = []\n",
    "masks_future = []\n",
    "masks_fr = []\n",
    "\n",
    "for i in tqdm(range(len(ths_op) - 1)):\n",
    "    \n",
    "    if i < len(ths_op)  - 2:\n",
    "        mask_temp_curr_pos = (ths_op[i] <= data[op1]) & (data[op1] < ths_op[i+1])\n",
    "        mask_temp_fut_pos = (ths_op[i] <= data[op2]) & (data[op2] < ths_op[i+1])\n",
    "        mask_temp_fr_pos = (ths_op[i] <= data[fr_op]) & (data[fr_op] < ths_op[i+1])\n",
    "    else:\n",
    "        mask_temp_curr_pos = (ths_op[i] <= data[op1]) & (data[op1] <= ths_op[i+1])\n",
    "        mask_temp_fut_pos = (ths_op[i] <= data[op2]) & (data[op2] <= ths_op[i+1])\n",
    "        mask_temp_fr_pos = (ths_op[i] <= data[fr_op]) & (data[fr_op] <= ths_op[i+1])\n",
    "   \n",
    "    masks.append(mask_temp_curr_pos)\n",
    "    masks_future.append(mask_temp_fut_pos)\n",
    "    masks_fr.append(mask_temp_fr_pos)\n",
    "    \n",
    "    \n",
    "control_mask = abs(data[fr_op] - data[fr_op2]) < 0.05  #avoid noisy opinion changes\n",
    "\n",
    "all_res = []\n",
    "all_res_l = []  #for confidence intervals (left edge)\n",
    "all_res_r = []  #for confidence intervals (right edge)\n",
    "\n",
    "for i in tqdm(range(len(masks))):\n",
    "    \n",
    "    new_res = np.zeros((len(masks_future), len(masks_fr)))\n",
    "    new_res_aux_l = np.zeros((len(masks_future), len(masks_fr)))\n",
    "    new_res_aux_r = np.zeros((len(masks_future), len(masks_fr)))   \n",
    "    \n",
    "    for j in range(len(masks_future)):\n",
    "        \n",
    "        for k in range(len(masks_fr)):\n",
    "            \n",
    "            \n",
    "            nom = sum(masks[i] & masks_future[j] & masks_fr[k] & control_mask)\n",
    "                \n",
    "            denom = sum(masks[i] & masks_fr[k] & control_mask)  \n",
    "            \n",
    "            \n",
    "            \n",
    "            if denom>0:\n",
    "                p = nom / denom\n",
    "                (conf_l, conf_r) = proportion_confint(nom, denom, alpha=0.05, method='wilson')\n",
    "            else:\n",
    "                p = None\n",
    "                (conf_l, conf_r) = (None, None)\n",
    "\n",
    "            new_res[j, k] = p\n",
    "            new_res_aux_l[j, k] = conf_l\n",
    "            new_res_aux_r[j, k] = conf_r  \n",
    "            \n",
    "    \n",
    "    all_res.append(new_res)\n",
    "    all_res_l.append(new_res_aux_l)\n",
    "    all_res_r.append(new_res_aux_r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "assigned-surge",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save slices\n",
    "#Note: len(ths_op)-1 = m\n",
    "files = os.listdir()\n",
    "    if 'Transition matrices' not in files:\n",
    "        os.mkdir('Transition matrices')\n",
    "\n",
    "for i in range(len(all_res)):\n",
    "    np.save(f'Transition matrices/Matrix_{len(ths_op)-1}_{i+1}.npy', all_res[i])  #if save transition matrix drawn from first iteration\n",
    "    #np.save(f'Transition matrices/Matrix_{len(ths_op)-1}_{i+1}_second_iteration.npy', all_res[i])  #if save transition matrix drawn from second iteration\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "excessive-china",
   "metadata": {},
   "source": [
    "# Display transition matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "based-diving",
   "metadata": {},
   "source": [
    "# $m = 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "former-treasury",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([[0.975, 0.025],\n",
       "        [0.952, 0.048]]),\n",
       " array([[0.066, 0.934],\n",
       "        [0.049, 0.951]])]"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = 2\n",
    "all_res = []\n",
    "for i in range(m):\n",
    "    all_res.append(np.load(f'Transition matrices/Matrix_{m}_{i+1}.npy'))  #first iteration\n",
    "    #all_res.append(np.load(f'Transition matrices/Matrix_{m}_{i+1}_second_iteration.npy'))  #second iteration\n",
    "P = [np.round(all_res[i].T, 3) for i in range(m)]\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sealed-strategy",
   "metadata": {},
   "source": [
    "# $m = 3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "hearing-senegal",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([[0.96 , 0.04 , 0.   ],\n",
       "        [0.942, 0.057, 0.001],\n",
       "        [0.907, 0.091, 0.002]]),\n",
       " array([[0.039, 0.952, 0.008],\n",
       "        [0.021, 0.969, 0.01 ],\n",
       "        [0.02 , 0.944, 0.036]]),\n",
       " array([[0.001, 0.082, 0.916],\n",
       "        [0.001, 0.07 , 0.929],\n",
       "        [0.001, 0.054, 0.945]])]"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = 3\n",
    "all_res = []\n",
    "for i in range(m):\n",
    "    all_res.append(np.load(f'Transition matrices/Matrix_{m}_{i+1}.npy'))  #first iteration\n",
    "    #all_res.append(np.load(f'Transition matrices/Matrix_{m}_{i+1}_second_iteration.npy'))  #second iteration\n",
    "P = [np.round(all_res[i].T, 3) for i in range(m)]\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "northern-hayes",
   "metadata": {},
   "source": [
    "# $m = 10$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "yellow-google",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([[0.942, 0.038, 0.003, 0.009, 0.009, 0.   , 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.938, 0.044, 0.005, 0.006, 0.005, 0.   , 0.003, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.945, 0.042, 0.007, 0.003, 0.002, 0.   , 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.947, 0.037, 0.008, 0.004, 0.003, 0.001, 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.939, 0.036, 0.009, 0.006, 0.005, 0.003, 0.001, 0.001, 0.   ,\n",
       "         0.   ],\n",
       "        [0.924, 0.043, 0.009, 0.008, 0.007, 0.004, 0.002, 0.001, 0.001,\n",
       "         0.   ],\n",
       "        [0.925, 0.045, 0.009, 0.006, 0.008, 0.002, 0.002, 0.002, 0.   ,\n",
       "         0.   ],\n",
       "        [0.908, 0.065, 0.006, 0.003, 0.   , 0.009, 0.003, 0.006, 0.   ,\n",
       "         0.   ],\n",
       "        [0.896, 0.067, 0.037, 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.881, 0.048, 0.   , 0.   , 0.   , 0.   , 0.048, 0.024, 0.   ,\n",
       "         0.   ]]),\n",
       " array([[0.13 , 0.766, 0.069, 0.017, 0.017, 0.   , 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.121, 0.776, 0.078, 0.016, 0.007, 0.002, 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.096, 0.821, 0.053, 0.022, 0.006, 0.001, 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.082, 0.839, 0.055, 0.016, 0.006, 0.001, 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.071, 0.839, 0.061, 0.019, 0.009, 0.002, 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.07 , 0.82 , 0.072, 0.023, 0.011, 0.002, 0.001, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.089, 0.762, 0.096, 0.032, 0.019, 0.002, 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.103, 0.761, 0.09 , 0.019, 0.024, 0.   , 0.003, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.114, 0.754, 0.079, 0.035, 0.018, 0.   , 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.159, 0.75 , 0.068, 0.   , 0.023, 0.   , 0.   , 0.   , 0.   ,\n",
       "         0.   ]]),\n",
       " array([[0.02 , 0.123, 0.733, 0.087, 0.037, 0.   , 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.011, 0.094, 0.785, 0.066, 0.037, 0.004, 0.002, 0.   , 0.   ,\n",
       "         0.002],\n",
       "        [0.015, 0.102, 0.792, 0.065, 0.022, 0.004, 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.009, 0.087, 0.823, 0.055, 0.022, 0.003, 0.   , 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.008, 0.075, 0.822, 0.065, 0.024, 0.005, 0.001, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.007, 0.073, 0.796, 0.082, 0.031, 0.009, 0.001, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.012, 0.088, 0.752, 0.1  , 0.034, 0.011, 0.001, 0.002, 0.   ,\n",
       "         0.   ],\n",
       "        [0.015, 0.096, 0.739, 0.097, 0.041, 0.01 , 0.001, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.025, 0.095, 0.714, 0.126, 0.025, 0.01 , 0.005, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.081, 0.757, 0.108, 0.054, 0.   , 0.   , 0.   , 0.   ,\n",
       "         0.   ]]),\n",
       " array([[0.002, 0.015, 0.099, 0.791, 0.074, 0.015, 0.003, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.006, 0.015, 0.087, 0.804, 0.073, 0.015, 0.   , 0.001, 0.   ,\n",
       "         0.   ],\n",
       "        [0.002, 0.018, 0.085, 0.811, 0.075, 0.009, 0.001, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.002, 0.011, 0.075, 0.838, 0.066, 0.006, 0.001, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.002, 0.008, 0.065, 0.837, 0.076, 0.011, 0.001, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.002, 0.007, 0.071, 0.809, 0.093, 0.015, 0.002, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.002, 0.009, 0.076, 0.788, 0.099, 0.021, 0.003, 0.001, 0.001,\n",
       "         0.   ],\n",
       "        [0.001, 0.015, 0.102, 0.758, 0.094, 0.021, 0.007, 0.001, 0.001,\n",
       "         0.   ],\n",
       "        [0.002, 0.012, 0.103, 0.763, 0.099, 0.017, 0.004, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.007, 0.042, 0.12 , 0.732, 0.092, 0.007, 0.   , 0.   , 0.   ,\n",
       "         0.   ]]),\n",
       " array([[0.002, 0.001, 0.011, 0.08 , 0.854, 0.051, 0.002, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.001, 0.004, 0.009, 0.067, 0.873, 0.04 , 0.004, 0.001, 0.001,\n",
       "         0.   ],\n",
       "        [0.001, 0.004, 0.01 , 0.07 , 0.874, 0.038, 0.003, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.001, 0.002, 0.008, 0.06 , 0.896, 0.031, 0.001, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.001, 0.001, 0.005, 0.054, 0.895, 0.041, 0.002, 0.001, 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.001, 0.005, 0.059, 0.871, 0.059, 0.004, 0.   , 0.   ,\n",
       "         0.   ],\n",
       "        [0.001, 0.001, 0.005, 0.068, 0.837, 0.078, 0.008, 0.001, 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.002, 0.008, 0.068, 0.827, 0.083, 0.011, 0.003, 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.001, 0.005, 0.066, 0.826, 0.09 , 0.009, 0.003, 0.   ,\n",
       "         0.001],\n",
       "        [0.   , 0.003, 0.013, 0.09 , 0.79 , 0.092, 0.005, 0.005, 0.   ,\n",
       "         0.003]]),\n",
       " array([[0.   , 0.   , 0.004, 0.005, 0.119, 0.816, 0.05 , 0.006, 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.   , 0.005, 0.01 , 0.119, 0.815, 0.045, 0.006, 0.   ,\n",
       "         0.   ],\n",
       "        [0.001, 0.001, 0.002, 0.012, 0.105, 0.838, 0.038, 0.004, 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.001, 0.003, 0.01 , 0.106, 0.85 , 0.029, 0.001, 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.   , 0.002, 0.006, 0.093, 0.866, 0.03 , 0.002, 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.   , 0.001, 0.006, 0.094, 0.848, 0.046, 0.003, 0.   ,\n",
       "         0.   ],\n",
       "        [0.   , 0.001, 0.002, 0.006, 0.095, 0.822, 0.07 , 0.005, 0.001,\n",
       "         0.   ],\n",
       "        [0.001, 0.001, 0.002, 0.007, 0.095, 0.803, 0.082, 0.008, 0.001,\n",
       "         0.   ],\n",
       "        [0.001, 0.001, 0.002, 0.009, 0.105, 0.789, 0.084, 0.009, 0.001,\n",
       "         0.001],\n",
       "        [0.   , 0.   , 0.007, 0.009, 0.093, 0.759, 0.119, 0.012, 0.   ,\n",
       "         0.   ]]),\n",
       " array([[0.002, 0.   , 0.   , 0.002, 0.009, 0.122, 0.8  , 0.062, 0.002,\n",
       "         0.   ],\n",
       "        [0.   , 0.002, 0.   , 0.002, 0.018, 0.131, 0.793, 0.05 , 0.004,\n",
       "         0.   ],\n",
       "        [0.001, 0.001, 0.001, 0.005, 0.016, 0.118, 0.804, 0.05 , 0.003,\n",
       "         0.003],\n",
       "        [0.001, 0.   , 0.001, 0.005, 0.013, 0.115, 0.822, 0.041, 0.002,\n",
       "         0.   ],\n",
       "        [0.   , 0.   , 0.001, 0.002, 0.01 , 0.104, 0.842, 0.038, 0.002,\n",
       "         0.   ],\n",
       "        [0.   , 0.   , 0.001, 0.002, 0.009, 0.1  , 0.836, 0.05 , 0.002,\n",
       "         0.   ],\n",
       "        [0.   , 0.   , 0.001, 0.002, 0.01 , 0.101, 0.809, 0.074, 0.004,\n",
       "         0.   ],\n",
       "        [0.001, 0.   , 0.001, 0.002, 0.012, 0.101, 0.783, 0.094, 0.005,\n",
       "         0.001],\n",
       "        [0.   , 0.001, 0.001, 0.005, 0.015, 0.101, 0.783, 0.086, 0.007,\n",
       "         0.001],\n",
       "        [0.003, 0.   , 0.   , 0.   , 0.021, 0.103, 0.782, 0.088, 0.003,\n",
       "         0.   ]]),\n",
       " array([[0.   , 0.   , 0.   , 0.   , 0.004, 0.022, 0.146, 0.774, 0.053,\n",
       "         0.   ],\n",
       "        [0.   , 0.   , 0.   , 0.003, 0.007, 0.014, 0.111, 0.815, 0.049,\n",
       "         0.   ],\n",
       "        [0.001, 0.   , 0.001, 0.003, 0.006, 0.014, 0.113, 0.824, 0.038,\n",
       "         0.001],\n",
       "        [0.   , 0.   , 0.001, 0.002, 0.005, 0.015, 0.112, 0.824, 0.038,\n",
       "         0.003],\n",
       "        [0.   , 0.   , 0.   , 0.001, 0.003, 0.011, 0.099, 0.85 , 0.034,\n",
       "         0.001],\n",
       "        [0.   , 0.   , 0.   , 0.001, 0.003, 0.01 , 0.092, 0.848, 0.044,\n",
       "         0.001],\n",
       "        [0.   , 0.   , 0.001, 0.001, 0.004, 0.011, 0.092, 0.83 , 0.06 ,\n",
       "         0.001],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.005, 0.01 , 0.091, 0.81 , 0.081,\n",
       "         0.002],\n",
       "        [0.   , 0.   , 0.   , 0.001, 0.004, 0.011, 0.104, 0.796, 0.085,\n",
       "         0.   ],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.003, 0.016, 0.101, 0.804, 0.075,\n",
       "         0.   ]]),\n",
       " array([[0.   , 0.   , 0.   , 0.   , 0.   , 0.012, 0.029, 0.118, 0.818,\n",
       "         0.024],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.005, 0.005, 0.022, 0.077, 0.863,\n",
       "         0.027],\n",
       "        [0.   , 0.   , 0.   , 0.005, 0.   , 0.008, 0.01 , 0.098, 0.844,\n",
       "         0.035],\n",
       "        [0.   , 0.   , 0.001, 0.003, 0.002, 0.006, 0.01 , 0.098, 0.851,\n",
       "         0.03 ],\n",
       "        [0.001, 0.   , 0.001, 0.001, 0.003, 0.004, 0.008, 0.087, 0.871,\n",
       "         0.023],\n",
       "        [0.   , 0.   , 0.   , 0.001, 0.001, 0.003, 0.007, 0.082, 0.88 ,\n",
       "         0.026],\n",
       "        [0.   , 0.   , 0.   , 0.001, 0.002, 0.004, 0.008, 0.082, 0.87 ,\n",
       "         0.034],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.003, 0.006, 0.008, 0.076, 0.871,\n",
       "         0.035],\n",
       "        [0.001, 0.001, 0.   , 0.   , 0.001, 0.004, 0.011, 0.095, 0.827,\n",
       "         0.06 ],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.   , 0.008, 0.004, 0.054, 0.888,\n",
       "         0.046]]),\n",
       " array([[0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.018, 0.   , 0.125,\n",
       "         0.857],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.   , 0.018, 0.018, 0.   , 0.161,\n",
       "         0.804],\n",
       "        [0.   , 0.   , 0.005, 0.   , 0.005, 0.005, 0.   , 0.014, 0.096,\n",
       "         0.876],\n",
       "        [0.001, 0.   , 0.   , 0.   , 0.006, 0.005, 0.005, 0.005, 0.081,\n",
       "         0.897],\n",
       "        [0.   , 0.   , 0.   , 0.001, 0.002, 0.004, 0.004, 0.005, 0.064,\n",
       "         0.919],\n",
       "        [0.001, 0.   , 0.   , 0.001, 0.001, 0.002, 0.005, 0.005, 0.076,\n",
       "         0.909],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.002, 0.002, 0.006, 0.007, 0.079,\n",
       "         0.904],\n",
       "        [0.001, 0.   , 0.001, 0.   , 0.   , 0.003, 0.006, 0.005, 0.074,\n",
       "         0.91 ],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.005, 0.   , 0.002, 0.014, 0.079,\n",
       "         0.9  ],\n",
       "        [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.018, 0.071,\n",
       "         0.912]])]"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = 10\n",
    "all_res = []\n",
    "for i in range(m):\n",
    "    all_res.append(np.load(f'Transition matrices/Matrix_{m}_{i+1}.npy'))  #first iteration\n",
    "    #all_res.append(np.load(f'Transition matrices/Matrix_{m}_{i+1}_second_iteration.npy'))  #second iteration\n",
    "P = [np.round(all_res[i].T, 3) for i in range(m)]\n",
    "P"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
